import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 常数定义
H_if_cm = 1  # 单位：cm^-1
h = 6.62607015e-34  # 普朗克常数，单位：J·s
c = 2.998e10  # 光速，单位：cm/s
H_if = H_if_cm * h * c  # 转换为J

T = 298  # 温度，单位：开尔文
h_bar = 1.0545718e-34  # 约化普朗克常数，单位：J·秒
k_B = 1.380649e-23  # 玻尔兹曼常数，单位：J/K

eV_to_J = 1.60218e-19  # eV到J的转换因子

# 创建ΔG和λM的网格，单位为eV
ΔG_range_eV = np.linspace(0, 0.3, 200)  # ΔG范围，从-2到2 eV
λM_range_eV = np.linspace(0.001, 1, 150)  # λM范围，从0.01到2 eV
ΔG_eV, λM_eV = np.meshgrid(ΔG_range_eV, λM_range_eV)

# 转换为J
ΔG = ΔG_eV * eV_to_J
λM = λM_eV * eV_to_J

# 根据给定公式计算kT
prefactor = (np.pi / h_bar) * np.abs(H_if)**2
exponential = np.exp(-( (ΔG + λM)**2 ) / (4 * λM * k_B * T))
kT = prefactor * (1 / np.sqrt(np.pi * λM * k_B * T)) * exponential

# 绘图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(ΔG_eV, λM_eV, kT, cmap='viridis')

ax.set_xlabel('ΔEst (eV)')
ax.set_ylabel('λM (eV)')
ax.set_zlabel('krisc (s^-1)')

plt.title('3D Surface Plot of krisc as a Function of ΔEst and λM')
plt.show()
